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1  Introduction 

Vast  verity  of  geophysical  flows  occur  in  density  stratified  fluids  in  a  presence  of  sheared 
mean  velocity  profiles.  Under  such  conditions  it’s  been  observed  that  vigorous  turbulence 
leads  to  generation  of  mixed  layers  within  the  fluid,  e.g.  the  regions  of  homogeneous  density 
separated  from  each  other  by  strong  density  gradients  (interfaces).  These  layers  could  be  of 
large  vertical  scales  as  well  as  of  fine  scales.  An  example  of  such  flows  is  the  destratification 
of  the  upper  layer  of  the  ocean  under  the  action  of  the  surface  wind  stress  and  the  formation 
of  multiple  mixed  layers  which  are  observed  within  the  seasonal  thermocline  in  the  ocean 
as  well  as  in  fresh  lakes  [10].  Such  structures  in  a  turbulent  flows  are  thought  to  play 
a  crucial  role  in  the  determining  the  turbulent  fluxes  of  momentum,  buoyancy  and  other 
tracers  which  in  turn  affect  the  further  dynamics  of  the  flow.  Nonetheless,  the  problem  of 
turbulent  transport  remains  one  the  most  difficult  ones  in  geophysics. 

Previous  studies  investigated  mechanisms  acting  in  stratified  turbulent  flows  -  a  review 
of  those  could  be  found  in  [4],  The  formation  of  small  scale  mixed  layers  in  stratified  fluids 
driven  by  external  mechanical  sources  of  turbulence  were  investigated  by  theoretically  by 
Balmforth  et  al  [1],  and  experimentally  by  Park  et  al  [6].  In  the  present  study  the 
dynamics  of  layers  is  explored  within  a  turbulent  fluid  in  a  presence  of  a  sheared  mean 
flow.  A  set  of  experiments  were  performed  in  which  a  linearly  stratified  fluid  is  driven 
by  the  stress  applied  by  a  rotating  horizontal  disk  at  the  surface.  This  setting  leads  to 
the  formation  of  the  surface  mixed  layer  extending  from  the  the  disk  into  the  interior  of 
the  fluid.  Its  evolution  with  time  and  its  dependence  of  the  stratification  and  the  rotation 
of  the  disk  are  obtained.  Furthermore,  small  scale  mixed  layers  were  observed  in  below 
the  surface  mixed  layer.  They  have  a  complicated  dynamics  with  multiple  regimes,  which 
were  attempted  to  characterize.  The  description  of  the  experiment  and  the  behavior  of 
the  observed  mixed  layers  is  presented  in  section  (2).  Further,  a  phenomenological  model  is 
developed  that  explains  the  mechanism  under  which  the  secondary  mixed  layers  are  formed. 
It  parametrized  the  action  of  the  two  mechanisms:  shear  generation  of  the  turbulence  due 
to  Kelvin-Helmholtz  type  of  instabilities,  and  the  mean  flow  generated  turbulence  though 
the  vortex  scraping  mechanism.  The  resulting  gradient-type  mixing  model  predicts  the 
formation  of  the  instabilities  in  a  particular  parameter  regime  that  lead  to  mixed  layers 
formation.  The  instabilities  are  related  to  those  described  by  Phillips  [7]  and  Posmentier 
[8].  The  model  description,  its  analysis  and  numerical  simulations  are  presented  in  section 
(3).  The  summary  of  the  work  is  in  section  4. 
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Figure  1:  Experimental  set-up:  cylindrical  tank  filled  with  stratified  fluid,  rotating  disk  at 
the  surface  of  the  fluid,  conductivity  probe 

2  Experiments 

A  series  of  experiments  were  made  where  a  stratified  fluid  was  forced  by  a  rotating  horizontal 
disk  at  the  surface  with  the  objective  to  investigate  the  dynamics  of  the  formed  mixed  layers 
under  various  cases  of  stratification  and  the  disk  rotating  speeds.  The  experimental  setting 
is  similar  to  that  used  by  Davies  et  al  [3],  who  investigated  the  flow  generated  during  the 
spin  up  process  and  further  established  secondary  circulation. 

2.1  Experiment  arrangement 

Experimental  set-up  is  shown  on  Figure  1.  A  cylindric  tank  of  30  cm  height  and  30  cm 
diameter  was  filled  with  a  salt  stratified  water.  Only  linear  stratifications  were  considered, 
which  were  obtained  with  a  standard  double  bucket  technique  [5].  The  range  of  density 
variations  within  the  fluid  was  from  2%  to  20%.  A  horizontal  disk  with  the  diameter  of 
24  cm  is  located  at  the  surface  of  the  fluid,  with  its  axis  of  rotation  aligned  with  the  axis 
of  the  cylinder.  The  rotation  speed  of  the  disk  is  controlled  by  the  motor  (attached  above 
the  disk)  that  gives  a  range  of  0.5  -  10  rad/s.  Measurements  of  vertical  density  profiles 
were  obtained  using  a  conductivity  probe  which  was  taking  100  measurements  per  1  cm  of 
fluid,  moving  vertically  with  a  speed  of  3  mm/s.  Measurements  were  taken  only  when  the 
probe  was  moving  down  to  avoid  fluxes  of  water  along  the  probe  that  would  contaminate 
the  data.  The  probe  Reynolds  number  is  of  the  order  of  100  (  this  is  an  upper  bound,  which 
is  reached  at  the  edge  of  a  disk);  thus,  measurements  are  not  considerably  affected  by  the 
presence  of  the  probe  on  a  scales  considered  in  the  current  study  -  0{cm).  The  profiled 
depth  of  the  fluid  was  20  cm,  whereas  the  height  of  it  is  27 cm.  Density  was  calculated  based 
on  conductivity  with  a  calibration  using  a  3rd  order  polynomial  as  a  mapping  function.  To 
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Figure  2:  Images  of  the  typical  interfaces  between  the  well  mixed  layer  and  a  linearly 
stratified  layer  below  it. 


reduce  temperature  effects  on  conductivity  the  fluid  allowed  to  reach  the  room  temperature 
(20°C)  before  the  start  of  the  experiment.  Thus,  the  control  parameters  that  define  the 
outcome  of  the  experiment  were  the  stratification  of  the  fluid  (expressed  as  a  Brunt -Vaisala 
frequency  N)  and  the  rotating  speed  of  the  disk  fi. 

There  are  multiple  non-dimensional  parameters  involved  in  current  experiments.  The 
most  relevant  one  is  the  Richardson  number,  which  will  be  defined  here  as  Ri  =  .  This 

definition  should  not  be  confused  with  the  local  Richardson  number  yf  o  which  varies  with 
time  and  the  location  in  the  fluid.  The  explored  range  of  Ri  G  (0.15  —  2.3).  Another 
parameter  would  be  the  disk  edge  Reynolds  number:  Re  =  ^  ~  105,  where  U  is  the 
azimuthal  velocity  at  the  edge  of  the  disk,  L  is  the  radius  of  the  disk,  and  /r  is  the  kinematic 
viscosity  of  water.  A  flow  with  such  high  values  of  Reynolds  number  could  be  considered 
turbulent  and  independent  from  the  exact  value  of  Re.  The  Schmidt  number  for  the  salt 
stratified  water  is  a  =  -  ~  700,  where  k  is  the  salt  diffusivity.  The  tank  aspect  ratio  as  well 
as  the  ration  of  the  disk  radius  to  the  height  of  the  water  column  is  ~  1.  All  geometrical 
parameters  as  well  as  the  Schmidt  number  were  fixed  throughout  all  experiments.  Thus,  the 
only  non-dimensional  parameter  that  describes  the  different  outcomes  of  the  experiments 
is  the  Richardson  number. 

2.2  Observations  of  density  profiles 

After  the  tank  is  filled  with  the  stratified  fluid  the  disk  is  set  in  motion.  As  was  mentioned 
previously,  the  disk  Reynolds  number  is  big  enough  that  the  flow  adjacent  to  it  becomes 
turbulent.  Such  flows  are  very  effective  in  diffusing  buoyancy  and  momentum.  Thus,  a 
mixed  layer  forms  below  the  disk,  with  almost  homogeneous  density  and  velocity  distribu¬ 
tions.  The  growth  of  such  a  layer  however  is  suppressed  by  the  stable  stratification  of  the 
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density  profiles 


Figure  3:  Density  profiles  plotted  for  different  times  after  the  start  of  the  experiment. 
Profiles  are  shifted  along  the  X-coordinate  by  0.03 g/cm5  from  each  other. 


fluid.  As  the  mixed  layer  deepens  the  gravitational  potential  energy  of  the  fluid  increases 
because  mixing  process  effectively  lifts  the  center  of  the  mass  of  a  water  column  by  mixing 
dense  waters  upwards  against  the  gravitational  field.  Furthermore,  mixing  in  a  stratified 
fluid  leads  to  formation  of  a  density  jump  at  the  base  of  the  mixed  layer.  The  density 
jump  suppresses  turbulence  and  effectively  separates  the  turbulent  flow  in  the  upper  mixed 
layer  from  almost  quiescent  layer  underneath  it.  Within  the  mixed  layer  there  is  a  mean 
azimuthal  flow  due  to  the  coherent  rotation  of  the  disk.  This  flow  is  also  suppresses  when  it 
reaches  strong  stratification  at  the  base  of  the  layer  creating  a  strong  shear  and  providing 
other  mechanisms  of  the  generation  of  turbulence  which  would  be  discussed  further. 

Figure  2  shows  images  of  two  interfaces.  The  top  one  appears  at  early  stages  in  the 
development  of  the  mixed  layer  and  is  characterized  by  weak  density  jump  across  it,  ener¬ 
getic  eddies  and  large  thickness  (few  cm).  The  one  shown  at  the  bottom,  forms  during  later 
stages  and  has  low  energy  eddies  of  smaller  sizes  (less  then  a  cm)  and  much  stronger  density 
jump  across  it.  The  signature  of  these  interfaces  is  clearly  present  in  the  vertical  density 
profiles,  the  time  evolution  of  which  is  shown  on  Figure  3.  Initially,  linearly  stratified  fluid 
get  mixed  from  the  top  and  forms  a  mixed  layer  as  well  as  the  density  jump  at  the  base 
of  it,  both  of  which  are  growing  in  magnitude.  The  interface  separating  the  mixed  layer 
from  almost  undisturbed  bottom  layer  changes  its  characteristics  as  well.  The  profile  taken 
at  t  =  0.2 h  (hours)  shows  the  presence  of  unstable  stratification,  which  resulted  from  the 
overturning  of  energetic  eddies  at  the  base  of  the  mixed  layer.  The  density  interface  for  this 
profile  corresponds  to  a  type  shown  on  Figure  2  (top).  At  the  profile  taken  at  t  =  6 h  it  is 
hard  to  see  the  density  structures  resulted  from  overturning  of  eddies  because  turbulence  is 
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Exp# 

N 

n 

Ri 

a 

A 

9 

0.5 

1.13 

0.2 

0.2 

0.42 

13 

0.5 

1.34 

0.14 

0.18 

0.5 

14 

0.52 

0.72 

0.52 

0.2 

0.32 

15 

1.37 

1.7 

0.65 

0.23 

0.27 

17 

1.3 

0.95 

1.87 

0.246 

0.19 

18 

1.1 

1.5 

0.54 

0.22 

0.29 

19 

1.12 

1.7 

0.43 

0.2 

0.36 

20 

1.18 

1.13 

1.09 

0.25 

0.21 

21 

1.3 

1.34 

0.95 

0.24 

0.23 

22 

0.98 

0.95 

1.05 

0.24 

0.24 

23 

0.84 

2.08 

0.16 

0.21 

0.46 

24 

0.78 

0.72 

1.18 

0.24 

0.23 

25 

1.46 

1.13 

1.67 

0.25 

0.18 

26 

1.7 

1.7 

1 

0.225 

0.22 

27 

1.7 

1.13 

2.25 

0.225 

0.24 

28 

1.72 

3.14 

0.3 

0.18 

0.41 

29 

1.6 

1.94 

0.67 

0.22 

0.29 

30 

1.42 

1.13 

1.57 

0.23 

0.22 

31 

1.5 

1.18 

1.6 

0.24 

0.22 

32 

0.94 

1.5 

0.4 

0.23 

0.29 

Table  1:  The  table  shows  the  experiment  settings:  initial  stratification  of  the  fluid  (N),  the 
disk  angular  velocity  Q.  The  last  two  columns  show  the  coefficients  of  the  power  law  fit  for 
evolution  of  the  mixed  layer  depth  h  =  Ata.  The  values  are  organized  and  sorted  by  the 
experiment  number. 


highly  suppressed  there  and  eddies  are  of  smaller  scale;  this  interface  corresponds  the  one 
on  Figure  2(bottom). 

The  presence  of  the  density  jump  and  a  velocity  shear  across  it  points  to  a  possibility 
of  generation  of  the  Kelvin-Helmholtz  type  of  instabilities  at  the  interface  that  leads  to  its 
further  erosion.  However,  these  type  of  instabilities  would  be  present  if  the  local  Richardson 
number  (IV2 /Uz2)  is  sufficiently  small.  Eventually  the  flow  would  reach  conditions  stable  to 
K-H  instabilities  as  the  velocity  of  the  mixed  layer  would  gradually  slow  down,  whilst  the 
density  jump  across  the  interface  would  increase.  At  this  point,  the  erosion  of  the  interface 
would  continue  due  to  the  mean  flow  acting  frictionally  on  turbulent  eddies  at  the  interface 
(vortex  scraping).  Thus,  whether  it’s  the  mean  flow  or  its  shear,  there  is  a  mechanism  which 
would  generate  mixing  at  the  stable  interface  and  lead  to  its  further  deepening. 

2.3  Development  of  the  primary  mixed  layer 

As  was  described  earlier,  a  mixed  layer  forms  at  the  disk  and  expands  into  the  fluid  reaching 
the  thickness  of  the  order  of  the  tank  size.  It  will  be  referred  to  as  the  primary  mixed  layer. 
In  certain  cases,  there  are  other  layers  that  form  below  it  (secondary  mixed  layers)  which 
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Figure  4:  A  non-dimensional  depth  of  the  primary  mixed  layer  is  plotted  as  a  function  of 
non-dimensional  time  (left)  and  on  a  log- log  plot  (right). 


are  of  smaller  scales  and  have  much  weaker  interfaces.  These  will  be  discussed  in  the  next 
section. 

For  many  geophysical  applications  it  is  important  to  know  the  growth  of  the  primary 
mixed  layer  with  time  as  it  controls  the  distribution  of  tracers  as  well  as  the  dynamics  of 
the  flow  .  To  present  the  observations,  the  mixed  layer  depth  h  is  non-dimensionalized  by 
the  depth  of  the  fluid  ( h  =  h/Hw),  and  time  t  by  the  period  of  rotation  of  the  disk  (t  =  Qt). 
Figure  4(left)  shows  the  time  evolution  of  the  primary  mixed  layer  for  a  typical  experiment: 
the  layer  grows  rapidly  at  first,  then  the  growth  slows  down  due  to  the  increasing  density 
jump  at  the  base  of  the  layer.  Figure  4(right)  shows  the  thickness  of  the  layer  as  a  function 
of  time  on  a  log-log  plot  with  a  corresponding  linear  fit  ( R 2  =  0.89).  This  fit  suggests  a 
power  law  h  =  Ata,  where  A  and  a  presumably  are  functions  of  the  Richardson  number  - 
the  control  parameter  of  the  experiment.  Indeed,  Figure  5  shows  that  there  is  a  particular 
dependence  of  the  fit  coefficients  on  the  Richardson  number.  The  power  a  stays  relatively 
constant  at  a  value  of  0.22  with  some  trend  towards  smaller  values  at  lower  Richardson 
numbers.  The  coefficient  A  of  the  fit  seems  to  scale  with  the  Richardson  number  as  A  ~ 
Ri~ 0.33  ag  ghown  in  Figure  5(right).  Thus,  the  experimental  data  shows  that  the  non- 
dimensional  mixed  layer  depth  scales  as  h  ~  i??'-0'33  f022. 

To  understand  the  obtained  dependence  of  the  mixed  layer  depth  evolution,  consider 
the  energy  balance  of  the  system.  The  source  for  the  generation  of  the  mean  flow  as  well 
as  the  mixing  of  the  stratified  fluid  is  a  power  produced  by  frictional  stress  due  to  motion 
of  the  disk.  It  scales  as  P  ~  tU  ~  U3,  where  r  ~  U2  is  the  frictional  stress  and  U  is  the 
characteristic  velocity  of  the  mixed  layer  near  the  disk.  Further,  a  linearly  stratified  fluid 
that  was  mixed  to  a  depth  h  would  have  its  gravitational  potential  energy  increased  by 
AGPE  =  J°h(p-  p0)gzdz  r\_/  N2h3,  where  po  is  the  initial  linear  density  profile  and  p  is  the 
one  obtained  by  adiabatic  mixing  the  initial  profile  within  the  layer  of  thickness  h.  At  last, 
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N2«l2  loci  (N2/rf) 


Figure  5:  Coefficients  of  the  power  fit  for  the  depth  of  the  mixed  layer  are  plotted  as 
functions  of  Ri.  On  the  left  is  the  coefficient  a  =  a(Ri );  coefficient  A  =  A(Ri )  plotted  on 
a  log- log  scale  (right). 


the  kinetic  energy  of  the  mixed  layer  would  scale  as  KE  ~  U2h  and  the  dissipation  D  is 
assumed  to  be  proportional  to  the  energy  production  ( D  =  7 P).  Thus,  the  energy  balance 
is  described  by  the  following  equation: 

Y(KE  +  PE )  =  P~D =7  ^-( U2h  +  N2h3)  ~  U3 

At  this  point,  an  assumption  is  made  that  the  flow  reaches  a  state  where  the  kinetic 
energy  of  the  mixed  layer  approaches  a  constant  value  and  the  energy  input  redistributes 
between  potential  energy  and  dissipation.  If  KE  =  cont  then  U  ~  /i-1/2  and  the  energy 
balance  model  becomes  h2ht  rsj  h  3/2  =7-  h  ^2/9  cphe  obtained  power  law  for  the  growth  of 
the  mixed  layer  is  consistent  with  observations  since  they  indicate  a  power  of  0.22  (Figure  5). 
The  dependence  of  the  fitted  coefficient  A  ~  RE1/3  suggests  that  the  assumed  constant 
value  of  the  kinetic  energy  is  a  function  of  Ri  as  well.  This  relation,  could  not  be  determined 
based  on  the  simple  energy  balance  model.  The  assumptions  made  in  constructing  an  energy 
balance  could  be  verified  with  measurements  of  the  velocity  of  the  fluid  using,  for  example, 
PIV  techniques.  In  current  experiment  these  measurements  were  not  produced,  but  for 
further  progress  there  are  of  vital  importance. 

2.4  Dynamics  of  a  secondary  mixed  layer 

Layers  of  weak  stratification  that  form  below  the  primary  layer  will  be  referred  to  as  the 
secondary  mixed  layers.  Figure  6  shows  an  example  of  a  developed  secondary  mixed  layer 
-  it  has  a  distinct  signature  in  density  profile  as  well  as  in  Brunt-Vaisala  frequency,  where 
bases  of  layers  (interfaces)  have  strong  stratification  and  the  layer  itself  is  weekly  stratified. 
Thus,  on  Figure  6(right)  the  strong  peak1  corresponds  to  the  base  of  a  primary  mixed  layer, 
and  the  peak  with  smaller  magnitude  corresponds  to  the  base  of  a  secondary  mixed  layer. 

1  Buoyancy  frequency  was  obtained  from  smoothed  density  profile  to  remove  noise  -  instrumental  as  well 

as  due  to  small  scale  overturning.  This  leads  to  broader  peaks  with  smaller  magnitudes;  nonetheless,  a  clear 

signature  of  mixed  layers  is  preserved. 
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Figure  6:  A  typical  density  profile  containing  a  primary  mixed  layer  as  well  as  a  secondary 
one  (left). Brunt-Vaisala  frequency  showing  peaks  at  the  bases  of  mixed  layers. 


Note  that,  secondary  mixed  layer  is  only  partially  mixed  as  apposed  to  the  primary  one 
that  is  completely  destratified.  Thus,  secondary  layers  are  characterized  by  the  magnitude 
of  a  peak  in  buoyancy  frequency  at  the  base  of  the  layer,  and  by  a  layer  thickness,  which  is 
the  distance  between  the  two  peaks. 

As  opposed  to  primary  mixed  layer,  that  monotonically  grows  with  time,  the  secondary 
layer  has  a  very  rich  dynamics.  Figure  7  shows  how  the  buoyancy  frequency  evolves  with 
time.  Here,  there  is  a  clear  indication  of  both  mixed  layers  -  the  primary  one  continuously 
grows  and  the  secondary  one  forms  and  decays  multiple  times  throughout  the  experiment. 
Characteristic  thickness  of  the  secondary  layer  is  of  few  centimeters.  The  evolution  process 
could  be  decomposed  in  several  stages.  In  the  beginning,  there  is  a  strong  turbulent  flow  that 
drives  the  growth  of  the  primary  layer  and  does  not  allow  the  formation  and  persistence  of 
the  fine  scale  secondary  layer.  Then  follows  a  period  of  time  when  the  secondary  layer  forms 
and  decays  with  a  characteristic  time  of  few  tens  of  minutes.  The  layer  forms  spontaneously 
having  a  finite  thickness  which  then  decreases  until  this  layer  merges  with  the  primary  one 
(collision  instability)  and  the  process  repeats.  The  next  stage  is  described  by  a  change  in 
the  way  the  layer  forms.  They  appear  at  the  base  of  a  primary  mixed  layer  having  small 
thickness  and  start  to  grow  building  up  the  density  jump  at  the  interface.  However,  the 
magnitude  of  a  jump  vanishes  with  time  (erosion  instability),  the  layer  disappears  and  the 
next  one  forms.  Finally,  comes  the  stage  when  the  layers  either  no  longer  form  or  stay  in  a 
locked  position  relative  to  the  primary  interface. 

The  evolution  of  the  secondary  mixed  layer  is  highly  dependent  on  the  Richardson 
number  for  the  experiment  and  in  a  current  time.  In  some  cases,  certain  stages  of  the 
evolution  could  be  skipped  which  makes  it  hard  to  classify  all  the  regimes.  Nonetheless,  an 
attempt  to  plot  a  regime  diagram  was  made.  Figure  8  shows  the  stages  of  the  evolution 
of  the  mixed  layer  for  different  experiments  as  a  function  of  non  dimensional  time  and 
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Figure  7:  Examples  of  the  time  evolution  of  the  square  of  the  Brunt-Vaisala  frequency  for 
different  values  of  the  Richardson  number  ( Ri  =  0.54, 1,  2.25, 1.67). 


Figure  8:  Regime  diagram  showing  different  stages  of  evolution  of  the  secondary  mixed 
layer 
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Richardson  number.  Below  the  red  curve  steps  generate  below  the  mixed  layer  and  decay 
by  merging,  whereas  in  between  the  red  and  blue  curves  mixed  layers  generate  from  the 
primary  mixed  layer  and  decay  by  erosion.  Above  the  blue  line  the  staircases  are  either  in 
locked  state  or  do  not  exist.  The  diagram  does  not  show  any  specific  dependence,  however 
this  does  not  imply  that  there  is  no  dependence  whatsoever.  There  is  an  obvious  lack  of 
data:  a  total  of  20  experiments  were  made  and  for  less  then  half  of  them  the  secondary  mixed 
layers  could  be  clearly  classified.  Furthermore,  since  there  are  no  theoretical  predictions  it 
is  not  clear  how  to  present  such  a  diagram:  the  most  obvious  choice  of  presenting  in  terms 
of  Richardson  number  and  time  might  not  be  the  clear  one  to  look  at. 

The  dynamics  of  the  secondary  mixed  layers  is  extremely  complicated.  Understanding  it 
to  a  level  of  making  predictions  requires  a  precise  theoretical  model  bucked  up  by  extensive 
observations  of  not  only  density  but  also  velocity  fields,  turbulent  buoyancy  and  momen¬ 
tum  fluxes.  However,  due  to  time  constrains  these  elaborate  measurements  could  not  be 
performed.  Thus,  the  theoretical  part  of  this  research  would  be  focused  on  explaining  the 
formation  of  secondary  layers  in  an  idealized  phenomenological  framework. 

3  Phenomenological  model 

Secondary  mixed  layers  arise  in  a  stratified  fluid  that  is  subjected  to  the  action  of  a  turbulent 
flow  with  a  mean  shear.  The  turbulent  nature  of  the  flow  would  be  the  key  to  understanding 
the  phenomenon.  The  proposed  mechanism  is  as  follows.  The  sheared  mean  flow  driven 
by  the  rotating  disk  is  a  source  of  turbulent  motion  which  inhomogeneously  distributes 
within  a  stratified  fluid  depending  on  its  stratification.  The  buoyancy  itself  is  affected  by 
turbulent  fluxes.  Thus,  the  formation  on  the  layers  could  be  related  to  the  instability  of 
the  turbulent  flux-buoyancy  relation  as  was  first  proposed  by  Phillips  [7]  and  Posmentier 
[8].  Suppose,  the  evolution  of  buoyancy  b  is  controlled  by  the  divergence  of  the  downward 
turbulent  buoyancy  flux  F: 

db  _  dF 

dt  dz 

and  this  flux  is  a  function  of  the  buoyancy  gradient  F  =  F(bz).  The  equation  has  obvious 
steady  state  solutions  with  constant  flux,  however  such  solutions  would  be  unstable  if  the 
flux  would  be  a  decreasing  function  of  buoyancy  gradient  e.i.  <  0.  Figure  9  shows  the 
mechanisms  of  the  instability.  A  linear  steady  state  buoyancy  profile  (dashed  line)  has  a 
value  of  buoyancy  gradient  that  corresponds  to  the  negative  derivative  of  the  buoyancy  flux. 
After  introducing  a  small  perturbation  (black  line)  there  will  appear  anomalous  buoyancy 
fluxes  (blue  arrows).  In  regions  where  the  perturbations  have  higher  buoyancy  gradient 
the  anomalous  flux  is  negative  and  in  it  is  positive  in  regions  of  lower  gradients.  As  could 
be  seen  from  the  evolution  equation  the  current  distribution  of  the  buoyancy  fluxes  would 
lead  to  a  growth  of  the  initial  perturbation,  meaning  that  the  base  state  is  unstable.  This 
instability  would  lead  to  buoyancy  profiles  consisting  of  a  series  of  mixed  layers  separated 
by  strong  interfaces  (staircases).  In  the  case  of  >  0  the  created  fluxes  would  lead  to  the 
decay  of  the  buoyancy  anomaly  -  this  state  is  stable. 

This  model  considers  the  buoyancy  equation  only  and  assumes  a  particular  shape  of 
the  flux-gradient  curve.  However,  the  exact  form  of  this  curve  would  depend  on  particular 
mixing  mechanisms  that  create  the  turbulent  buoyancy  fluxes.  In  the  following  section 
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Figure  9:  Schematic  representation  of  the  instability  mechanism  that  leads  to  layer  forma¬ 
tion.  Dashed  line  represents  the  initial  buoyancy  profile,  black  line  is  the  slightly  perturbed 
profile,  anomalous  buoyancy  fluxes  denoted  by  blue  arrows. 

a  gradient-type  mixing  model  would  be  devised  that  would  describe  the  development  of 
buoyancy  profiles  under  the  action  of  the  mean  flow  generated  turbulence.  The  formation 
of  the  layers  in  the  model  would  be  related  to  the  mechanism  described  above  by  calculating 
the  exact  form  of  the  flux-gradient  relation.  The  model  would  consist  of  evolution  equations 
written  for  buoyancy  as  well  as  for  turbulent  kinetic  energy  that  is  drained  from  the  mean 
flow.  It  is  aimed  to  describe  the  secondary  mixed  layers  only,  so  it  is  assumed  that  the 
primary  layer  has  been  formulated  and  the  mean  velocity  profile  is  established. 

3.1  Model  formulation 

The  approach  taken  in  constructing  a  mixing  model  is  similar  to  that  implemented  by 
Balmforth  [1]  and  Barenblatt  [2].  The  model  parametrizes  turbulent  vertical  transports 
as  gradient  transports  using  eddy  diffusivity  and  viscosity.  Thus,  the  conservation  laws  for 
horizontally  averaged  buoyancy  b  and  turbulent  kinetic  energy  e  are  as  follows 


db 
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dt 

~dz 

'KWz 

de 

d 
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'KWz 
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where  S  is  a  source  term  for  the  energy  and  k,  is  the  turbulent  diffusivity.  The  source 
term  represents  production  and  destruction  of  the  turbulent  energy  and  in  this  formulation 
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would  consist  of  four  terms: 


1 


S  =  —nbz  —  a—  H - kUz2  H - — — U2e1/l2l 

k  1  +  Ri  1  +  Ri 


-l 


(3) 


As  it  was  mentioned  previously,  the  mixing  of  buoyancy  leads  to  an  increase  in  a  gravita¬ 
tional  potential  energy  of  the  system.  Thus,  the  role  of  the  first  term  (— k62)  is  to  convert  the 
eddy  kinetic  energy  into  potential  energy  preventing  the  artificial  accumulation  of  the  total 
energy  E  =  f^(e  —  bz)dz.  The  second  term  —  cc^-  is  the  dissipation,  which  is  on  dimensional 
grounds  proportional  to  the  energy  e  and  to  the  inverse  of  the  eddy  turnover  time  r  ~  Ac/e; 
a  is  a  small  parameter  which  controls  the  intensity  of  the  dissipation  and  effectively  defines 
the  number  of  turnover  periods  over  which  an  eddy  would  dissipate.  The  last  two  terms 
are  the  production  terms  due  to  the  action  of  a  mean  flow.  From  visual  observations  of  the 
experiments  there  seem  to  be  two  mechanism  in  action.  At  low  local  Richardson  numbers 
(Ri  =  N2/Uz2),  where  the  destabilizing  effect  of  the  shear  is  dominating  the  stabilizing 
effect  of  the  stratification,  the  turbulence  is  generated  due  to  Kelvin-Helmholtz  type  insta¬ 
bilities.  These  instabilities,  however,  are  suppressed  at  high  Ri.  Hence  the  representation  of 
the  energy  production  is  of  the  form  ^^nU2.  At  high  Ri  the  turbulence  is  generated  by 
the  mean  flow  acting  on  turbulent  eddies  (vortex  scraping)  by  inducing  a  frictional  stress. 
The  amount  of  power  introduced  by  such  a  mechanism  is  proportional  to  U2el^2l~1 ,  where 
U2  is  scaling  for  the  stress  exerted  by  the  mean  flow,  e1/2  is  a  characteristic  velocity  of  an 
eddy  and  l  is  its  length  scale.  Factors  1^j  and  are  chosen  to  ensure  that  the  two 
mechanisms  would  act  in  their  corresponding  regime  of  Richardson  numbers. 

Eddy  diffusivity,  as  well  as  viscosity,  is2  proportional  to  le1^2.  To  close  the  system  the 
characteristic  eddy  size  l  needs  to  be  defined.  Since  the  model  aims  to  describe  mixing  in  the 
stratified  fluid  there  are  several  length  scales  involved.  In  the  limit  of  strong  stratification 
the  size  of  eddies  would  be  proportional  to  the  Ozmidov’s  length  scale  (e/bz)1/2.  This 
restriction  comes  from  the  fact  that  eddies  overturn  density  layers  leading  to  increase  in 
potential  energy.  Thus,  for  a  given  energy  of  an  eddy  and  the  ambient  stratification  there 
would  be  a  restriction  on  its  size.  In  the  limit  of  weak  stratification,  the  characteristic 
length  scale  for  an  eddy  l  =  d  should  be  related  to  the  geometrical  parameters  e.i.  the  size 
of  the  tank3.  Thus,  to  interpolate  the  two  length  scales  in  the  limits  of  strong  and  weak 
stratifications  expression  for  the  eddy  length  scale  l  is  chosen  to  be 


l 2  d2  '  e/bz 


(4) 


where  7  is  a  parameter  that  determines  the  length  scale  in  the  limit  of  strong  stratification. 
It  should  be  pointed  out  that  the  model  is  written  in  terms  of  locally  defined  dynamical 
variables,  meaning  that  all  the  scales  are  controlled  by  the  local  dynamics  and  do  not  depend 
on  global  profiles  of  the  variables.  However,  this  doesn’t  always  take  place  in  reality.  For 
example,  the  model  used  for  the  eddy  length  scale  states  that  in  the  case  of  locally  weak 
stratification  the  length  should  be  equal  to  d  -  some  geometrical  restriction.  However, 
vertical  profiles  of  stratification  containing  mixed  layers  bounded  by  strong  stratification 

2  For  simplicity  eddy  viscosity  and  diffusivity  are  taken  to  be  the  same,  however  this  assumption  could 
be  relaxed 

3A  restriction  of  the  length  scale  in  weakly  stratified  case  avoids  singularities  in  the  model. 
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would  impose  an  equivalent  ’geometrical’  constrain  on  an  eddy  size  in  the  location  within 
the  mixed  layer  -  the  eddy  size  could  not  be  greater  than  the  thickness  of  a  mixed  layer  with 
strong  interfaces.  Thus,  the  parameter  d  should  be  considered  as  effective  one  that  takes 
into  account  geometrical  constrains  as  well  as  the  global  buoyancy  profile.  In  the  current 
model,  for  simplicity  it  is  assumed  to  be  constant. 

A  resulting  model  is  effectively  a  set  of  coupled  nonlinear  diffusion  equations,  which 
require  boundary  and  initial  conditions  to  be  completed.  As  boundary  conditions,  the 
fluxes  of  energy  and  buoyancy  are  set  to  be  zero  at  the  top  and  the  bottom  boundaries  of 
the  domain  e.i.  ez  =  0,  bz  =  0  at  z  =  0,  H.  The  initial  conditions  would  be  discussed  in  a 
section  3.3.  Thus,  the  complete  equation  set  of  the  mixing  model  is  the  following: 
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where  a,  7,  d  are  fixed  parameters;  characteristic  velocity  U  and  velocity  shear  Uz  rep¬ 
resent  the  forcing  for  this  system  and  assumed  to  be  constant.  In  principle,  equations  6 
to  9  could  be  combined  to  give  an  evolution  equation  for  the  turbulent  diffusivity  k  -  the 
variable  that  determines  the  evolution  of  the  buoyancy  profile. 

Since  the  purpose  of  the  model  is  to  describe  the  formation  of  mixed  layers  it  is  conve¬ 
nient  to  use  the  buoyancy  gradient  (g  =  bz)  instead  of  buoyancy  as  a  dynamical  variable. 
The  mixed  layers  would  have  a  distinct  signature  in  buoyancy  gradient  profile:  low  val¬ 
ues  within  the  mixed  layer  and  sharp  peaks  at  the  interfaces.  Further,  it  is  useful  to  use 
non-dimensional  equation  set  with  the  following  definition  of  dimensionless  quantities: 
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In  dimensionless  form  the  equations  are 
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b.c. 

:  ez  =  0,  g  =  0  at  z  =  0,  H 

(14) 

where  a  parameter  Q  describes  the  relative  importance  of  the  mean  flow  and  the  mean 
shear  in  the  energy  production  term;  the  local  Richardson  number  is  proportional  to  the 
buoyancy  gradient  only  (Ri  =  g/ 7). 
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Figure  10:  Steady  state  buoyancy  flux  plotted  as  a  function  of  the  buoyancy  gradient  for 
different  values  of  the  forcing  parameter  Q. 


3.2  Steady  states  and  their  stability 

Consider  a  steady  state  of  the  model,  which  has  a  constant  energy  and  buoyancy  gradient 
(eo>5o)  =  const.  It  violates  the  boundary  conditions,  nonetheless  away  from  the  influence 
of  the  boundaries  it  is  useful  to  consider  its  stability.  Such  steady  states  should  have  no 
production  of  turbulent  kinetic  energy  (S  =  0)  in  order  for  dynamic  variables  not  to  change 
with  time: 


e20  ,  7«o  +  Qg0 

-n0go  -  ay - 1 - — - f—  =  0 

« o  1  +  50/7 


(15) 

(16) 


This  sets  a  constrain,  which  effectively  defines  energy  and  buoyancy  flux  as  a  function 
of  buoyancy  gradient  (or  Richardson  number)  for  the  considered  steady  state  solutions: 


e0  =  E(g0)  (17) 

«o5o  =  F(go)  (18) 

As  was  discussed  earlier,  the  instabilities  would  develop  if  the  initial  profile  has  the 

Richardson  number  which  is  within  the  range  of  the  negative  slope  of  the  buoyancy  flux. 

Figure  10  shows  the  steady  state  buoyancy  flux  as  a  function  of  Richardson  number  (Ri  = 
g/ 7).  The  numerical  values  of  parameters  are  chosen  as  a  =  0.1, 7  =  0.2.  The  representative 
mean  flow  parameter  Q  =  0.8,1.5,2.56.  Curves  with  values  of  Q  <  2.56  show  a  non¬ 
monotonic  behavior,  having  a  region  where  the  derivative  of  the  flux  is  negative.  At  low 
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values  of  Ri,  the  instabilities  are  not  possible  due  to  extremely  high  turbulent  energy  created 
by  shear  instabilities  that  are  not  suppressed  by  the  stratification.  At  high  Richardson 
numbers  the  buoyancy  flux  is  controlled  by  the  turbulence  produced  by  the  mean  flow 
which  is  proportional  to  g 4/2,  hence  the  increase  of  the  buoyancy  flux  with  the  increase  in 
stratification.  However,  the  shear  production  term  is  ~  e/53/2,  which  means  that  there  is 
strong  penalty  in  energy  production  for  an  increase  in  stratification  -  hence  the  turbulent 
flux  due  to  shear  production  should  decrease  with  the  stratification  at  high  Richardson 
numbers.  Thus,  the  effect  of  the  mean  shear  tends  to  destabilize  the  system  and  create 
mixed  layers,  whereas  the  mean  flow  turbulence  has  a  stabilizing  effect.  Hence,  in  the 
intermediate  range  of  Richardson  the  slope  has  a  negative  slope  due  to  shear  production 
term,  however  this  region  is  of  finite  size  due  to  the  action  of  the  mean  flow  turbulence  at 
higher  values  of  Ri.  This  means  that  the  instability  could  not  grow  unbounded  in  time:  at 
some  point  the  effect  of  the  mean  flow  turbulence  would  suppress  the  growth.  Thus,  the 
existence  and  development  of  instabilities  depend  crucially  on  the  relative  importance  of 
the  mean  flow  to  shear  production  terms  which  is  determined  by  the  value  of  the  parameter 
Q.  If  the  forcing  Q  is  greater  than  2.56  the  flux- gradient  curve  monotonically  increases  with 
the  value  of  Richardson  number  -  stabilizing  effect  of  the  mean  flow  is  dominant.  For  such 
forcing  all  the  buoyancy  profiles  would  be  stable. 

The  described  instabilities  could  now  be  identified  from  the  formulated  model  equations 
and  their  growth  rates  could  be  quantified.  Linearizing  the  equations  around  the  base  state 
(eo;5o)  and  assuming  exponentially  growing  solutions  for  perturbations  (e,g)  ~  exp(Xt  — 
ikz )  leads  to  a  linear  stability  problem: 
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where  f  =  Kg  is  the  buoyancy  flux,  S  is  the  source  term  in  the  energy  equation  and  all  the 
partial  derivatives  are  calculated  at  the  point  (eo,  go)-  Thus,  for  a  given  value  of  parameter  Q 
and  initial  value  of  the  Richardson  number  a  growth  rate  A  could  be  calculated  as  a  function 
of  the  wavenumber  k.  Figure  11  shows  growth  rates  calculated  for  Q  =  1.5  and  for  values 
of  base  state  Ri  =  0.3, 0.7, 1.3.  The  growth  rate  is  positive  for  Ri  =  0.7  which  is  located  at 
the  negatively  sloped  part  of  the  flux  curve  (Figure  10).  The  growth  rate  is  negative  for  the 
values  of  Richardson  numbers  of  0.3  and  1.3  that  belong  to  the  monotonically  increasing 
part  of  the  slope  -  these  are  stable  solutions.  An  important  property  of  the  written  model  is 
that  in  the  case  of  the  unstable  steady  states  the  growth  rate  of  the  perturbations  has  has 
a  high  wavenumber  cutoff4  which  means  that  small  scale  perturbations  do  not  grow  with 
time  -  instead,  they  decay.  This  property  allows  to  perform  direct  numerical  simulations 
of  the  development  of  the  instabilities  and  their  evolution  leading  to  the  formation  of  the 
mixed  layers. 


3.3  Numerical  simulations 

Numerical  simulations  of  the  mixing  model  were  performed  using  first  order  in  time  and 
second  order  in  space  numerical  scheme  with  an  adaptive  time  stepping.  Computational 
domain  is  z  =  [0, 1200]  having  600  uniformly  distributed  grid  points.  Initial  conditions  for 

4The  high  wavenumber  cutoff  is  due  to  the  diffusion  of  the  turbulent  kinetic  energy 
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Figure  11:  Growth  rate  plotted  as  a  function  of  the  wavenumber  for  different  values  of  the 
base  state  Richardson  numbers  and  for  a  fixed  parameter  Q  =  1.5  (as  predicted  by  the 
linear  stability  analysis). 


Figure  12:  Numerically  simulated  evolution  of  the  instabilities  manifested  in  the  vertical 
profiles  of  the  Richardson  number. 
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Figure  13:  A  snapshot  of  model  variables  at  t  =  le4:  Richardson  number  (left),  buoyancy 
flux  (middle)  and  turbulent  kinetic  energy  (right). 


the  numerical  simulation  are  chosen  to  be  Ri  =  0.8  for  the  buoyancy  profile5,  e  =  0.01  for 
the  turbulent  kinetic  energy  and  the  forcing  parameter  Q  =  1.  Based  on  a  linear  stability 
analysis  such  a  configuration  should  be  unstable.  Indeed,  the  numerical  simulations  confirm 
it.  Figure  12  shows  snapshots  of  the  modeled  buoyancy  gradient  at  different  time  -  it  is 
clearly  seen  that  the  instabilities  are  developing  in  the  form  of  sharp  peaks  in  buoyancy 
gradient  profiles.  First,  the  instabilities  develop  close  to  boundaries  because  there  are  larger 
perturbation  from  the  steady  state  their.  Then,  they  appear  in  the  interior  and  eventually 
fill  the  whole  domain  except  for  small  regions  adjacent  to  boundaries  which  are  affected  by 
zero  flux  boundary  conditions. 

Figure  13  shows  a  snapshot  of  the  buoyancy  gradient  g,  turbulent  kinetic  energy  e  and 
a  corresponding  buoyancy  flux  /.  The  buoyancy  flux  seems  to  be  relatively  constant  in 
the  interior  region,  however  having  isolated  regions  of  low  flux.  The  energy  distribution  is 
higher  in  the  regions  of  weak  stratification  and  has  minimums  at  within  highly  stratified 
interfaces.  The  instabilities  develop  having  a  wavelength  corresponding  to  the  one  that 
has  a  maximum  growth  rate  predicted  by  a  linear  stability  analysis.  However,  it  has  been 
observed  that  after  the  initial  development  stage  the  structure  tends  to  coarsen  with  time 
by  merging  or  decaying  mechanisms.  The  individual  peaks  are  quasi-steady  solutions  of  the 
equations  and  evolve  by  weakly  interacting  with  each  other.  The  result  of  the  interaction 
could  be  considered  as  a  secondary  instability  that  leads  to  the  coarsening  through  erosion 
or  collision  of  the  adjacent  interfaces  [9].  The  time  scales  for  the  two  instabilities  would 
depend  on  the  base. 

Figure  14  shows  the  time  evolution  of  the  buoyancy  gradient  and  the  buoyancy  flux 
throughout  the  whole  simulation.  The  initial  development  stage  lasts  till  the  time  t  ~  4000 
and  is  characterized  by  the  formation  of  the  peaks  in  buoyancy  gradient.  The  overall 
structure,  however,  tends  to  coarsen.  The  peaks  either  merge  with  each  other  creating  a 

5To  satisfy  initial  conditions  at  the  boundaries  the  distribution  of  Ri  was  chosen  such  that  it  is  constant 
within  the  interior  of  the  domain  and  continuously  approaches  zero  at  the  boundaries. 
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Figure  14:  Time  evolution  of  the  modeled  vertical  profiles  of  the  Richardson  number  (left) 
and  the  buoyancy  flux  (right). 


higher  and  broader  peaks,  or  they  decay  in  amplitude  and  disappear  amplifying  the  adjacent 
peaks(  Figure  14(top)  ).  The  flux  on  the  other  hand,  started  from  an  initial  condition  which 
has  a  weak  constant  value  within  the  interior  and  after  the  instabilities  have  developed  its 
value  also  remained  constant,  however  at  a  higher  value  (Figure  14(bottom)).  Note  that 
the  merging  and  decaying  events  are  accompanied  by  a  transient  decrease  in  buoyancy 
flux.  This  is  clearly  seen  at  Figurel3,  where,  for  example,  the  two  peaks  at  z  ~  550 
or  z  ~  750  are  about  to  merge  and  the  corresponding  values  of  the  buoyancy  flux  are 
decreased.  The  edge  regions  which  are  characterized  by  weak  stratification  and  buoyancy 
flux  are  continuously  growing  free  of  instabilities.  The  end  result  of  such  a  simulation  would 
be  the  full  homogenization  of  the  buoyancy,  since  system  is  forced  by  the  mean  flow  which 
is  a  continuous  source  of  turbulent  mixing. 

The  numerical  simulation,  showed  that  a  series  of  mixed  layers  could  form  in  a  region 
subjected  to  the  action  of  the  mean  flow  turbulence.  However,  the  simulation  was  performed 
in  an  idealized  setting  so  that  it  could  confirm  the  results  of  the  linear  analysis  of  the  stability 
of  the  model.  The  model  only  conceptually  shows  the  mechanism  of  layer  formation.  In 
reality,  there  are  many  more  factors  that  would  contribute  to  the  evolution  of  the  layers. 
First  of  all,  the  mixing  model  uses  constant  characteristic  shear  and  the  mean  flow  as  a 
mechanisms  for  the  generation  of  turbulence.  However,  the  mean  flow  is  affected  by  changes 
in  stratification  and  should  be  accounted  by  the  model  if  the  instabilities  develop  at  a  similar 
time  scale  as  the  variations  in  mean  flow.  Furthermore,  in  the  experiments  there  is  a 
nonuniform  buoyancy  gradient  profile,  as  it  has  a  dominant  peak  at  the  base  of  the  primary 
mixed  layer.  The  existence  of  this  layer  affects  the  evolution  of  the  adjacent  secondary 
mixed  layers.  Nonetheless,  the  model  shows  the  important  features  of  the  evolution  of  the 
secondary  mixed  layers:  spontaneous  generation  of  the  layers,  the  merging  of  the  layers  and 
their  decay  -  all  of  which  have  been  observed  in  the  experiments.  Even  though  the  proposed 
model  simulates  similar  behavior  to  the  observations,  the  assumption  of  the  model  should 
be  verified  experimentally  -  this  would  require  elaborate  measurements  of  the  mean  flow  as 
well  as  the  of  the  turbulent  buoyancy  and  momentum  fluxes. 
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4  Conclusions 


A  series  of  experiments  were  performed  on  a  linearly  stratified  fluid  subjected  to  the  surface 
frictional  stress  induced  by  a  rotating  disk.  The  control  parameters  for  the  experiments  were 
the  magnitude  of  the  stratification  ( N 2)  and  the  rotating  speed  of  the  disk  (O)  resulting 
in  a  single  non-dimensional  parameter  Ri  =  N'2 /O2.  Two  types  of  mixed  layers  have  been 
observed  -  the  primary  layer  extending  from  the  surface  into  the  interior  of  the  fluid,  and 
the  secondary  mixed  layer  of  much  smaller  scale  that  forms  underneath  it.  The  growth 
of  the  primary  layer  with  time  seems  to  obey  a  power  law  h  ~  .  This  result  is 

consistent  with  the  energy  balance  model  that  assumes  a  constant  kinetic  energy  withing 
the  layer. 

The  secondary  mixed  layers  show  a  complicated  behavior,  having  several  regimes  de¬ 
pending  on  time  and  the  Richardson  number  of  the  experiment.  A  physical  mechanism 
for  the  formation  of  the  secondary  mixed  layers  has  not  been  clearly  identified  from  the 
observations  due  to  the  lack  of  observational  data,  namely  the  turbulent  fluxes  and  the 
mean  flow.  Nonetheless,  the  hypothesis  has  been  put  forward  that  relates  the  formation 
of  these  layers  to  an  instability  of  the  flux-gradient  type  that  could  arise  in  stratified  tur¬ 
bulent  flows.  A  one  dimensional  mixing  model  has  been  devised  and  the  instabilities  were 
investigated  theoretically  as  well  as  through  numerical  simulations.  The  resulting  mixed 
layers  form  only  in  a  particular  regime  that  puts  constraints  on  the  magnitude  of  the  mean 
flow  with  respect  the  its  shear  and  the  Richardson  number  for  the  profile.  Furthermore, 
the  dynamics  predicts  the  destruction  of  the  mixed  layers  through  merging  and  decay  -  the 
precesses  observed  in  experiments. 

Based  on  the  mixing  model,  the  results  of  the  experiments  could  be  interpreted  as 
follows.  Supported  by  the  disk  generated  turbulence,  the  primary  mixed  layer  grows  creating 
a  large  density  gradient  at  the  base  of  it.  At  the  same  time,  the  mean  flow  is  established 
within  the  primary  layer.  The  flow  is  suppressed  as  it  reaches  high  stratification  and  thus 
the  shear  is  created  at  the  base  of  the  mixed  layer.  This  sets  a  particular  value  of  the  model 
parameter  Q,  which  reflects  the  relative  importance  of  the  mean  flow  to  the  mean  shear. 
The  magnitude  of  the  flow  as  well  as  of  its  shear  changes  with  time  as  the  mixed  layer  grows, 
meaning  that  Q  is  a  function  of  time  as  well  as  space.  The  formation  of  the  secondary  mixed 
layers  below  the  primary  is  observed  when  the  necessary  conditions  for  the  instabilities  are 
met:  the  parameter  Q  is  smaller  then  the  critical  value,  and  the  local  Richardson  numbers 
belong  to  the  unstable  part  of  the  flux-gradient  curve.  Such  a  constrain  identifies  a  confined 
region  in  a  fluid  where  the  instabilities  could  occur.  The  instabilities  within  this  region  lead 
to  the  development  of  the  layers,  however  the  magnitude  of  the  instabilities  is  suppressed 
by  the  stabilizing  effect  of  the  mean  flow  turbulence.  That  explains  why  secondary  mixed 
layers  are  not  completely  destratified.  After  the  spontaneous  development  of  the  layers, 
comes  the  stage  of  slow  evolution  when  the  layers  either  merge  with  the  primary  interface 
of  decay  by  erosion  -  consistent  with  the  behavior  of  the  observed  layers.  The  change  in  the 
regime  is  the  shift  from  the  dominating  collision  instability  to  the  erosion  instability  which 
is  associated  with  the  change  in  the  base  state. 
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